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Abstract 

Following the Non-Relativistic QCD approach we use a gauge invariant 
smearing method with factorization to measure the excitation energies for 
a heavy QQ system on a X 48 lattice at /3 = 6.2. The results come from 
averaging over an ensemble of 60 QCD configurations. In order to enhance 
the signal from each configuration we use wall sources for quark propaga- 
tors. The quark Hamiltonian contains only the simplest non-relativistic 
kinetic energy term. 

The results are listed for a range of bare quark masses. The mass 
splittings are insensitive to this variable though there are a slight trends 
with increasing quark mass. For an appropriate choice of UV cut-off 
(a-^ = 3.2Gev) the mass spectrum compares reasonably well with the 
experimental values for the spin-averaged energy gaps of the T system. 

We also present results for the DE and DT waves for the lowest bare 
quark mass. The results are consistent with degeneracy between the two 
types of D wave. This encourages the idea that even with our simple quark 
Hamiltonian the departure from rotational invariance is not great. 



1 Introduction 



In a previous paper |Jl| we studied heavy quark bound states appropriate to 
a description of the J/ip andT systems using the non-relativistic approach of 
Lepage (NRQCD) 0, H, ^, |^. We investigated the lowest bound states for S", 
P and D waves ignoring spin effects for the quarks using gauge configurations 
from the UKQCD collaboration on a 16^ x 48-lattice with a /?-value of 6.2. In the 
present paper we report results on the first excitations in the 5* and P-channels for 
this system (again without spin). We obtained these results using 60 quenched 
QCD configurations from the UKQCD collaboration on a 24'^ x 48 lattice at 
(3 = 6.2 . The new results are reasonably consistent with our previous ones but 
considerably more precise. 

Our results are based on the construction of a number of smeared and un- 
smeared operators that couple to the appropriate channels and the measurement 
of their cross correlators. The smeared operators are constructed in a gauge 
invariant manner. Using a simple subtraction procedure we show that the corre- 
lation functions do indeed have a multi-exponential structure. Our best estimates 
of the lowest states and the first excited states in both the S and P-channels of the 
QQ system are established by performing consistent correlated fits to the mea- 
sured operator correlators using appropriately factorizing two-exponential forms. 
Some three-exponential fits were attempted to test the range of applicability of 
the fits but did not lead to different conclusions. 

2 Quark Propagators 

The quark propagator in a given gauge field background is 



where the angle brackets indicate averaging over the quark degrees of freedom 
{ip{x)} and X = (x, t), ?/ = (y, 0) . When t = 0, G{x, y) = 5xy • 

The evolution for the (non-relativistic) quark propagator G{x,y) is 



where t denotes a unit step in the time direction and n is the order of the time- 
step update as discussed by Thacker and Davies P]. We set G{x, y) = for t < . 
The modified update is necessary for stability at certain values of the bare quark 
mass. In this paper we use n = 3 . The hamiltonian Hq is that appropriate to 
non-relativistic propagation ignoring spin of a quark of mass M, namely 



G{x,y) = {ij{x)ij\y)) 
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The covariant finite differences A+, A are given by their usual expressions (we 
have suppressed all colour and spin indices) 

AtG(a;, y) = U^{x)G{x + fi,y)- G{x, y) , (4) 
ATG{x, y) = G{x, y) - ul{x - ft)G{x - jl,y) . (5) 

3 Smeared Operators 

The ideal method for detecting excitated states in a given channel is to construct 
operators each of which couples only to one of the states. An alternative approach 
is to accept as a starting point a basis set of operators with quantum numbers 
appropriate to the channel of interest and to recognize that an intermediate state 
will couple to each of these operators in a unique way so that the exponential 
contribution associated with that state to the cross correlators of the basis op- 
erators will have a factorizing form. This is the approach we have adopted in 
dealing with the gauge invariant smeared operators we construct and use in our 
simulation. 

The operators we investigated for the S'-channel were in addition to the stan- 
dard point operator (we use x{^) to denote the anti-quark degrees of freedom) 

C>(o)(x) =x^(2;)V'(2;)+h.c. , (6) 

a set of operators of the form 

OM(x) = Y.x\x) {Mf{x)^{x + mjl) + M^{x)i^{x - mfi)) + h.c. , (7) 
A 

where the /i-sum is over space hke directions and the matrices M^(a;) and M'^{x) 
have the (appropriately ordered) product forms 

m— 1 m 

Mj:{x) = n U[.{x + vfi) and Mf{x) = IJ ^lix - lyfi) . (8) 

For the P-channel we use a family of operators of the form 

Of'\x) = x\x) (^MJ^{x)^{x + mfi) - Mf{x)tlj{x - mfi)) + h.c. , (9) 
The DE wave operators are 

Ogf ^™\a;) = x\x) {Mf{x)i,{x + mfi) + Mf{x)^{x - mfi) 

-M^{x)7jj{x + mv) - M^{x)il){x - mi>)) + h.c. , (10) 
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and the DT wave operators are 

O^I^'^\x) = x\x) (Af ^Af) + Af^A^T)) V^(x) + h.c. , (11) 

where 

a(™V(x) = MJ^{x)^{x + mix) - M'J\x)^{x - mfi) . (12) 
The correlation functions we measure are 

for 5- wave analysis, and 

FLit) = E mi^)ony)) ' (14) 

for the P-wave analysis. Here V = 24^, the spatial lattice volume. 
All of the above operators are of the form 

0(x) = x^(x)vl>(x) , (15) 

where for an appropriate set of (S'f/(3)-matrix) coefficients {Cxx'} 

^{x) = Y.C^^>^{x') , (16) 

x' 

and x' = (x', t) . A typical correlation function can be expressed as 

Fu{t) = ^Y.{^^TG^2{x,y)G^{x,y)) + c.c. , (17) 

^ xy 

where 

Guix,y) = {¥'\x)¥'^Hy)) , (18) 
where 1 and 2 indicate the two (possibly the same) smeared operators. Of course 
the Green's function Gi2{x,y) = J^x' C^l'G'^'^K^' yV) where Q^^\x',y) can be cal- 
culated with an appropriate change of initial condition by the same method as 
the original quark Green's function. 

4 Wall Source Method 

Because the computing overhead is considerable it is desirable to extract as much 
signal as possible from each pass through a gauge field configuration. To this end 
we modify the measured correlators F^^{t), F^^{t) and F^(t)as follows. Eq([T7|) 
leads to a method of computation for our typical correlation function Fi2{t) 
that requires the evaluation of the Green's functions for at least a representative 
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sample of y- values on the initial time slice if we wish to maximize the information 
to be extracted from each configuration. This is computationally onerous. An 
alternative procedure is the following. We replace the single y-summation in 
eq ([T7|) with a double sum thus 

Fi2{t) = ^ E (TrGi2(a;,2/)Gt(x,2/')) +C.C. , (19) 
x,y,y' 

where y' = (0,y') . We now rely on the gauge field averaging to eliminate the 
contributions from the gauge non-invariant off-diagonal terms in the above double 
(y, y')-sum leaving only the contribution from the gauge invariant diagonal terms 
for which y = y' • The disadvantage of the method is that the off-diagonal 
contributions provide noise even if they do average to zero. The advantage of 
the method is that we pick up all the diagonal terms in one pass since it is only 
necessary to compute the objects of the form 

g{x) = Y.Gix,y) , (20) 
y 

which satisfies the same equation as G{x,y) and the initial condition 

^?(0,x) = l . (21) 
Similar remarks apply to ^12(3^) = J2y Gi2{x,y) . We have then 

Fu{t) = ^j:(^Tguix)g\x)) , (22) 

In practice we do find that the method works well and does provide a good signal 
relatively economically. It is implicit in the discussion that no gauge fixing has 
been imposed on the ensemble of gauge fields. However because of the limitations 
of the data set the averaging procedure may not work perfectly. The different 
treatment of the two operators in the correlation function may mean that the 
symmetry F^^'^{t) = F^^'^(t) no longer holds. This does not destroy factor- 
ization and we allow for the asymmetry in fitting the data. This procedure is 
similar to the multiple origin approach first utilised by Kenway0 and Billoire et 
al.p], except we seed all sites on the initial timeslice. A quark wall source was 
also used by Gupta et al.^ although they fix to Coulomb gauge. 

5 Data Fitting 

As indicated above we try to fit the correlation function F^^{t) with the multi- 
exponential form 

^lW = E7(a)™7(a)ne-*'("'* • (23) 
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Our main results are obtained using a two exponential form requiring eight pa- 
rameters. This allows us to obtain estimates for the lowest 5* and P states 
together with the first excitations. Our statistical method involved a correlated 
least squares fit based on an estimate of the complete set of variances and cross 
correlators of all the fitted quantities. Since our results are based on 60 statis- 
tically independent gauge field configurations we felt that eight parameters was 
a reasonable number for the fitted form. We did perform three exponential fits 
in certain cases with twelve parameters. Where these seemed reliable they were 
consistent with the two exponential fits but with considerably less tight errorbars. 
The results we quote are from separate S and P-channels fits. We also carried 
out a combined S and P-channel fit but obtained results that were little different. 

The limitations of the data led us to confine ourselves to two operators per 
channel in any one fit. The precise form of the 2x2 matrix of correlators was 




At) 

(24) 

This is equivalent to the form in eq(p3|). The asymmetry in the factorized forms 



is represented by the departure of the parameters p^^^ and p'-^-* from unity. Note 
that we have parametrized the splitting AM between the two levels explicitly 
since this is the quantity of direct interest. 

The basis of the fitting procedure is the estimate of the correlation matrix of 
results. At any one time these comprised the two direct and two cross correlators 
for two operators evaluated on 48 time slices. The correlation matrix was there- 
fore of dimension 192 x 192 . Our data is extracted from 60 independent gauge 
configurations. The correlation matrix is therefore of rank r < 60 and therefore 
necessarily singular. In practice the effective rank of the correlation matrix is 
even less than this since beyond a certain point the eigenvalues become so small 
their estimation from the data is not reliable. The least squares fitting procedure 
and the associated error estimates require the use of the inverse of the correlation 
matrix. It is necessary and indeed correct to restrict the inversion of the matrix 
to an appropriate subspace that is spanned by eigenvectors with eigenvalues large 
enough for reliable estimation from the data. The dimension of the subspace is 
referred to as the Singular Value Decomposition (SVD) cut. 

In assessing the results of the fitting procedure we examined cases with a range 
of values of initial off-set and SVD cuts for different combinations of smeared 
operators. Our criterion for a choice of result was that the x^-value be acceptably 
near unity per degree of freedom and that the error was the best (usually the 
first) of a range of reasonably good and statistically consistent fits. 
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6 Explicit Diagonalization Scheme 



Before exhibiting the resuhs of the correlated fits we show directly the existence 
of a second exponential by means of a diagonalization method. Liischer and Wolff 
1^ have shown that the eigenvalues of the correlation matrix are of the form 



(l + 0(e-^*'w*)) (25) 



where AM(q) is the distance of state M(a) from other states. Thus we evaluate 
the eigenvalues of the correlation matrix using the appropriate Numerical Recipes 
routines |]lT|. Fig. 1 shows the result for Ma = 1.5 for the 5- wave combination 
Oq (x) and Of{x) - the ground state S-wave is suppressed revealing the existence 
of the exponential associated with the first excited state. Fig. 2 shows effective 
mass plots for the 15* and 25 states obtained from these graphs. Figs. 3 & 4 show 
similar results for the IP and 2P states. The results are reasonably consistent 
with those of the correlated fits discussed below which were used to produce the 
quoted numbers. 

In order to obtain reasonably smooth plots the effective mass was defined as 
M(t) = 0.25 * log (^^1^^ , (26) 

where 

A{t) = (F(t) + wF{t + 1) + w^F{t + 4) + w''F{t + 3)) /4 , (27) 
and w is chosen to render the terms in the sum of comparable size. 



7 Correlated Fits for S and P Waves 



The results of the correlated fits are shown in Table 1 . Also shown are the 
operators used to obtain the results, the per degree of freedom, the SVD-cuts 
and the off-sets at which a reasonable statistical stability set in. 



Mo a 


Chan. 


Ma 


AM a 


xVdof 


SVD-cut 


Off-set 


Operators m 


1.5 


S 


1.1152(8) 


0.198(24) 


25.6/18 


26 


17 


and 4 


P 


1.243(12) 


0.116(12) 


19.2/20 


28 


11 


2 and 6 


2.0 


s 


0.9788(7) 


0.194(12) 


11.53/12 


20 


12 


1 and 4 


p 


1.112(13) 


0.105(20) 


14.48/20 


28 


11 


2 and 6 


3.0 


s 


0.8287(7) 


0.179(17) 


21.68/18 


26 


15 


and 4 


p 


0.969(17) 


0.090(17) 


11.6/14 


22 


11 


2 and 6 



Table 1: The results for the ground states and first excited splits in the 5* and 
the P channels. 
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It is also interesting to compare the various mass splits with the spin-averaged 
values of the T-system. We use a conversion factor from lattice units to physical 
units of = 3.2 . The results are listed in Table 2 . 
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It is clear that the pattern of mass sphtting for the lower two bare masses 
(Mga = 1.5 & 2.0) is reasonably close to the actual splitting for the T except for 
the 25 level which appears to be too high but exhibits a downward trend as M^a 
increases. The other splits are less sensitive to changes in the bare mass. 



Mo a 


1S-2S 


IS -IP 


IS-2P 


1P-2P 


1.5 


0.198(24) 
0.634(77) Gev 


0.128(12) 
0.410(38) Gev 


0.244(17) 
0.781(54) Gev 


0.116(12) 
0.371(38) Gev 


2.0 


0.194(8) 
0.621(26) Gev 


0.133(13) 
0.426(42) Gev 


0.238(24) 
0.762(77) Gev 


0.105(20) 
0.336(64) Gev 


3.0 


0.179(17) 
0.573(54) Gev 


0.140(17) 
0.448(54) Gev 


0.230(24) 
0.736(77) Gev 


0.090(17) 
0.288(54) Gev 


Expt(T) 


0.563 Gev 


0.430 Gev 


0.795 Gev 


0.365 Gev 



Table 2: The results for various mass splits in lattice and physical units compared 
to the spin-averaged results for the T-system. The conversion factor is = 3.2 
Gev. 

The ratios of mass splits is independent of the choice for . For each 
bare quark mass these ratios are listed in Table 4 taking the central value of the 
15* — IP split as the base. 



Moa 


1S-2S 


IS -IP 


1S-2P 


1P-2P 


1.5 


1.55(19) 


1.00(12) 


1.91(13) 


0.91(9) 


2.0 


1.46(6) 


1.00(13) 


1.79(18) 


.79(15) 


3.0 


1.28(12) 


1.00(17) 


1.64(17) 


0.64(12) 


Expt(T) 


1.31 


1.00 


1.85 


0.85 



Table 3: The results for various mass splits expressed as ratios to the central 
value of the IS* — IP split and compared to the spin-averaged results for the 
T-system. 

These results are not dissimilar to those of the phenomenological non-relativistic 



quark models |jT2|. The spectrum is relatively independent of the bare quark mass 
though the 25* state is rather high. Apart from the anomalously high 25" state 
the ratios that fit best correspond to a bare quark mass somewhere between 
Mqq = 1.5 and M^a = 2.0 . Using the same conversion factor as above we find 
these correspond to Mq = 4.8 Gev and Mq = 6.4 Gev. This is to be compared 
with the P-quark mass of ~ 5 Gev suggested by the mass of the T itself. If 
we take the IS* — IP mass split for Mqa = 1.5 as the correct basis on which to 
calculate then we find = 3.4(3) Gev which encompasses the above value. 
Our results for a^^ are higher than suggested by a measurement of the string 



tension a . At P = 6.2, aa^ = 0.026(1) Q, |T3| that is ^/aa = 0.161(3) . If we use 
the phenomenological value ^/a = .42 Gev we obtain = 2.6(1) Gev. This is in 
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line with other estimates of a~^. Another way of expressing this discrepancy is to 
note that our lattice calculation at /3 = 6.2 yields a ratio ^/o^/ AM{1S—1P) ~ 1.25 
whereas the phenomenological result is ~ 1.0 . The question then is whether or 
not there is a reasonable explanation of this discrepancy for our model. One 
answer is to recognize that the string tension is associated with the long range 
part of the quark potential while the 15" — IP split comes about as a result of a 
balance between long and short range effects in the potential. The short range 
force is controlled by the strong coupling evaluated at a higher momentum, g*, 
than that, q, associated with the string tension. The main difference between 
the quenched and unquenched theories is the differential renormalization of the 
strong coupling a{q) at a given q due to the vacuum polarization effects of light 
quarks. If we fix the string tension to be the same in both theories then we 
have au{q) = aqiq) However because the quenched coupling runs faster than the 
unquenched one we have au{q*) > otq^q*) . The short range force in the quenched 
theory will therefore be weaker than in the unquenched case. Because the (S-wave 
states are particularly sensitive to the short range part of the Q — Q force they will 
be more deeply bound in the unquenched theory than in the quenched one. The 
P-waves will be controlled more by the longer range part of the force associated 
with the string tension. This will tend to leave the P-waves unchanged between 
the two theories with the result that for a given string tension AM(15' — IP) will 
be greater in the unquenched relative to the quenched theory. In turn this will 
yield a lower value for the ratio \/a/ IS.M{1S — IP) for the unquenched relative 
to the quenched theory in line with our results. Similar results will hold for for 
any two quantities associated with different momentum scales. In the quenched 
theory they will yield different estimates for while unquenched theory (by 
definition) will produce consistent estimates. 



8 Correlated Fits for D Waves 

In Table 4 we show the results for a two exponential fit to the two versions of the 
D wave operators for the bare quark mass MqO = 1.5 . 



il7o« 


Clian. 


Ma 


AM a 


X-/(M 


SVD-cut 


Ofi-sot 


(Jperators m 


1.5 


DE 


1.373(11) 


0.315(13) 


23.7/18 


26 


6 


2 and 6 


DT 


1.388(19) 


0.389(41) 


22.5/14 


22 


6 


1 and 4 



Table 4: The results for the ground states and first excited splits in the DE and 
the DT channels. 



The encouraging feature of these results is the degeneracy within errors not 
only of the basic states in the two channels but also of the first excited states. 
The conclusion is that to a good approximation the cubical symmetry of the 
(spatial) lattice is replaced by rotatational symmetry. Expressed in physical units 
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AM{1S — ID) = 0.83(4) Gev if we again use a^^ = 3.2 Gev. There is so far no 
observed D wave for the T system but the corresponding state for charmonium 
{J/ip) has AM{1S — ID) = 0.702 Gev. Given the simple and approximate nature 
of our heavy quark Hamilton this is an encouraging result. The quality of the 
results from the simulation restricted the application of the fitting procedure 
to the range t < 20 so the outcome for the mass gap may be expected to be 
on the high side. This circumstance may also explain why the measured gap 
AM{1D — 2D) ~ 1 Gev is implausibly high. It is interesting that it showed in 
both versions of the D wave spectrum. 



9 Conclusions 

We have measured the QQ mass splittings for the radial excitations of the S 
and P waves using the non-relativistic heavy quark propagators calculated from 
quenched gluon configurations from the UKQCD collaboration. Our measure- 
ments were based on 60 QCD configurations on a 24^ x 48 lattice. The Hamil- 
tonian we used to calculate the quark propagators was of the simplest kind con- 
taining only the kinetic energy contribution and omitting all higher corrections. 
The results emerged from two-exponential correlated fits to pairs of smeared op- 
erator correlators. Our best results yielded statistical errors ~ 10% for the mass 
splittings. 

On a broad picture our results are not inconsistent with the pattern of spin- 
averaged splittings of the T-system. In particular the ratios AM{1S — IP) : 
AM{1S — 2P) : AM{1S — ID) seem roughly correct. The 2P wave shows a 
slight dependence on the bare quark mass. We have not yet determined the 
dependence of the ID wave on this mass. The absence of an experimental ID 
state for T restricts us to a comparison with the corresponding charmonium state 
or theoretical quark model predictions but that comparison is encouraging. If we 
base our evaluation of the cut-off strictly on the IS* — IP mass split then we find 
the value = 3.4 Gev. This very much in line with scaling predictions from 
results of the corresponding calculations performed at values of f3 =5.7 and 6.0 



1^ . This suggested a bare quark mass for the kinetic energy Hamiltonian in 
the range 4.8 to 6.4 Gev. 

However there are two obvious problems that present themselves. The first is 
that the ratio involving the string tension ^/a/ AM{1S — IP) is measured as 1.25 
compared to a phenomenological value of 1.0 . It is plausible that this discrepancy 
may be due to the overly strong running of the coupling in the quenched theory. 
The second is the anomalously high value of AM{1S — 2S)/AM{1S — IP) and its 
sensitivity to the bare quark mass. It may be that the deficiencies of the quenched 
approximation can resolve this problem also. An alternative explanation is that 
there is a measurement problem with the S wave channel. Future measurements 
using different smeared operators constructed in the Coulomb gauge from quark 
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model wave functions should help to resolve the issue. 
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Figure Captions 



Fig. 1 The results for the S'-wave propagator Fqq (t) (circles) together with 
the the results of the diagonalization procedure applied to the correlation 
matrix for the operators Oq^^ and 04"^ (diamonds) revealing the contribu- 
tion of the 25'-state. 

Fig. 2 The effective mass plot for the IS (upper graph) and 2S states (lower 
graph). The estimates obtained from the correlated fit are indicated by a 
full line. 

Fig. 3 The results for the P-wave propagator -Fi^H^) (circles) together 
with the the results of the diagonalization procedure applied to the cor- 
relation matrix for the operators 02^^ and Og^^ (diamonds) revealing the 
contribution of the 2P-state. 

Fig. 4 The effective mass plot for the IP (upper graph) and 2P states 
(lower graph). The estimates obtained from the correlated fit are indicated 
by a full line. 
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